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ABSTRACT 

We run mean-field shearing-box numerical simulations with a temperature- 
dependent resistivity and compare them to a reduced dynamical model. Our simula- 
tions reveal the co-existence of two quasi-steady states, a 'quiet' state and an 'active' 
turbulent state, confirming the predictions of the reduced model. The initial conditions 
determine on which state the simulation ultimately settles. The active state is strongly 
influenced by the geometry of the computational box and the thermal properties of 
the gas. Cubic domains support permanent channel flows, bar-shaped domains exhibit 
eruptive behaviour, and horizontal slabs give rise to infrequent channels. Meanwhile, 
longer cooling time scales lead to higher saturation amplitudes. 
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The observed accretion luminosity of astrophysical disks de- 
mands that angular momentum be transported at a rate far 
in ex cess of that possible by mole cular viscosity in a laminar 
flow (IShakura fc Svunvaev 19751). The m agneto-rotational 
instability (MRI) (|Balbus fc Hawlevlll99ll ) provides a more 
efficient mechanism: the MRI generates turbulent motions 
that are strongly correlated and which significantly increase 
the effective flux of angular momentum. It is hence regarded 
the most promising candidate responsible for the observed 
anomalous transport in magnetised disks. Our main obser- 
vational probes of disks exploit their radiation properties, 
which in turn strongly depend on their thermal and chemical 
properties. So in order to connect observations with the dy- 
namical state of accretion disks we must better understand 
the relationship between MRI turbulence and the chemistry 
and thermodynamics of accretion disks. The goal of this pa- 
per is to make progress to this end. 

We shall focus principally on the role of the thermo- 
dynamics and in particular, resistive dissipation. From the 
outset we assume a generic model of radiative loss and fast 
chemistry, so that the resistivity r\ depends only on temper- 
ature. The latter temperature dependence, however, gives 
rise to an interesting feedback on the turbulent dynamics. 
Fluctuations in temperature will not only alter the resis- 
tive scale, and hence the turbulent cascade of energy (see 
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IFromang et al J 120071 ; iLesur fc Longarettll2007l ). it will also 
influence the linear MRI modes which drive the turbulence 
itself in regions where the resistivity is sufficiently high. 

In order to clarify the intriguing relationship between 
the disk t hermodynamics and the MRI-induced turbulent 
dynamics, Balb us fc Lesa ffrc (2008) devised a simple 2 vari- 
able dynamical system, which describes the evolution of the 
temperature, on the one hand, and the magnitude of the 
turbulence, on the other. The model incorporates both the 
dynamics of the MRI, heating (coupling the turbulent fluc- 
tuations to the temperature), radiative cooling and a tem- 
perature dependent resistivity r\ (coupling the temperature 
to the turbulent fluctuations). Under very general circum- 
stances the simple system possesses two stable equilibria (or 
fixed points) . One fixed point corresponds to a cold flow — 
the 'quiet' state, in which the resistivity is high, the MRI 
shuts off, and the temperature is set by radiative balance. 
The other fixed point corresponds to a hot turbulent flow 
— the 'active' state, in which dissipative heating balances 
the radiative gains and losses. Which state the system se- 
lects ultimately depends on the initial conditions. Accretion 
disks are potentially susceptible to this kind of phase sep- 
aration. They need only experience a broad enough range 
of temperatures so that both active a nd quiet states coex - 
ist. For example, protop lanetary disks (IFromang et alj|2002l ) 
and dwarf-novae disks l|Gammie fc Menou 19981 ) may pos- 
sibly experience large enough resistivities to have a direct 
influence on outer-scale dynamics. 

In this paper we show that this framework can be used 
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to interpret more sophisticated models of disks, namely nu- 
merical simulations of shearing boxes. We restrict ourselves 
here to unstratified shearing boxes with a net mean verti- 
cal magnetic flux. We also investigate the effect of changing 
the aspect ratio of the computational box, as this seems to 
strongly influence the qualitative results. We first describe 
the reduced dynamical model (Section 2). The details of the 
numerical setup are presented in the Appendix. 

We then present our results, first comparing the simu- 
lations in cubic domains with the reduced model (Section 
3), then discussing the aspect ratio dependence of the sim- 
ulations' results (Section 4). In Section 5 we draw our con- 
clusions. 



2 SHEARING BOX REDUCED MODEL 

2.1 Governing equations for the 3D simulations 

Our simulations use a standard shearing-box setup. The 
frame of reference rotates at circular angular velocity f2. Ra- 
dial, azimuthal and vertical directions are labelled by local 
Cartesian coordinates x, y and z. The origin of the frame 
follows an unperturbed fluid element moving in a circular 
orbit. The radial logarithmic derivative of Q is 

1 dQ. 2 , 



* 2dln.fi 1 

and characterises the local shear. It takes the value q = 
3/2 for a Keplerian shear and we adopt this value for the 
remainder of the paper. This shearing sheet system is then 
supplemented with shearing box bou ndary conditions at the 
edges of the computational domain l|Hawlev et al.l[l995l ). 

The fundamental dynamical equations in this rotating 
frame are the mass continuity equation, 



|+v.(H = o, 



(i) 



where p is the mass density of the gas and v is its velocity, 
and the Navier-Stokes equation with a kinematic viscosity 
v, 
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where p is the thermal pressure, B is the magnetic field, 
$ = qQx 2 is the tidal potential and cry = \(piVj + djVi) — 
^dkVkSij is the stress tensor. Vertical gravity is neglected. 
The induction equation is 



|=Vx(«B-,(T)J) 



(3) 



where J = VxB and T](T) is the resistivity (a decreasing 
function of the temperature, see equation ()A3|) ). 

We run both isothermal and non-isothermal simula- 
tions. In the non-isothermal case, we adopt an ideal gas 
equation of state with adiabatic index 7 so that the internal 
energy is 



(4) 



and the equation for the Lagrangian derivative of the en- 
tropy reads 

e Dln ^ P 7) = nJ' 2 + pvv- V« - pA(T) - V-(XPVT) (5) 

where A(T) = aT + b is a linear net cooling function (in- 
cluding radiative gains and losses), \ is the thermal diffu- 
sion coefficient and T is the specific energy defined as e/p 
(it is hence proportional to the temperature). The adiabatic 
index 7 has been fixed to the standard value of 5/3 for the 
present study. More details on the thermal properties of our 
simulations are found in Appendix I A3I 

Equations Q to ((5J form the governing system of which 
we seek approximate numerical solutions. The steady state 
solution consists of the state of homogeneous density and 
linear shear: p — po and v = Axy = vq and we take a 
background magnetic field constant and vertical, B = Bo z. 
We denote the perturbation velocity and magnetic field as 
u — v — vo a nd b = B — Bp z. We use a version of the 
ZEUS3D code l|Stone fc Normanl [l99 2a,b) suited to our needs 
(Appendix lAl contains the full details of our numerical method). 



2.2 Simple dynamical model 

The reduced model of Balbus fc Lesaffre] l|2008t ) may be repre- 
sented as: 



^ = a(T)y - Ay n 
at 



(6) 



for the evolution of a generic turbulent amplitude y in the gas 
and 

dT 

_ = Wy 2 - A(T) (7) 
at 

for its temperature T. The time t is normalised to the orbital 
timescale 1/Q. The linear growth/damping rate of y is denoted 
by ct. It depend s on the temperatu re T via the resistivity rj(T). 
We follow iBalbus fc Lesaffre] d2008h and use 



a(T) = cr,[l- v (T)/r h ] 



(8) 



where cr t is the ideal MRI rate of growth at the computational box 
wavelength (for r\ = 0) and r;* is the critical resistivity for MRI 
at that wavelength (i.e. when r\ = r]„, the fastest growing mode in 
the computational box is marginally stable). The parameters A 
and n account for the non-linear saturation dynamics of the gas. 
The function A(T) denotes the radiative losses. And the constant 
W accounts for the dissipation of the turbulent kinetic energy at 
small scales. For stationary turbulence, this is equivalent to the 
energy input at large scales. 

But the comparison is quantitatively better if we improve 
equation (0 in our dynamical model to account for the finite 
time over which energy is transferred from large to small scales. 
We hence adopt the new system of reduced equations 

' b! ~a(T)y-Ay" (9) 



and 



at 



£ = w y-^-A(T). 

at at 



(10) 



The additional term dy 2 /dt with a time-derivative does not 
change the fixed points of the dynamical system and only slightly 
modifies their stability properties. 



2.3 Phase diagram of the reduced model 

IBalbus fc Lesaffre! j200Sh have shown that, under very general 
circumstances (namely: rj(T) decreasing and A(T) increasing to 
infinity), the above dynamical system possesses the same phase 
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diagram (illustrated on Fig. 1). It exhibits three fixed points. We 
have already discussed the two stable fixed points 'quiet' and 
'active' states in the introduction. We label them respectively 'Q' 
and 'A'. 

The third fixed point is alway s unstable and occurs for in- 
termediate temperatures: following Bal bus fc Lesaffrc (2008|) we 
label it T. An initial state (y,T) which finds itself near the co- 
ordinates of the point T ends up at 'A' or 'Q' depending on its 
position relative to the separatrices of the point T. These two 
separatrices define four quadrants in the (T, y) phase space. A 
cloud of various points evolving around point T splits into two 
cloudlets which end up one on the quiet region and one on the 
active region. In effect, this system describes the phase separation 
between active and quiet regions. 



2.4 Reduced variables 



Phase space 




„ 4 (to point Q) g . 

Temperature (T) 



In order to compare the numerical solution with the reduced 
model, we must extract from the simulations two quantities anal- 
ogous to T and y at each time step. We label these numerical 
variables T and y. 

We start with the total energy conservation (see equation (8) 
of lHawlev et al.lll995t for example) where we neglect the viscous 
stresses: 



±(e + m)=*-(W xy )-(pk{T)). 



(11) 



Here the angles denote volumic averages over the computational 
domain, e is the internal energy, m = \pv 2 + -g^b 2 + p$ is the 
total mechanical energy (including the tidal potential <3?) and 
W X y = pu x u y — j^b x b y is the total radial-azimuthal stress tensor 
of the perturbation. The variables p, v, u and b are respectively 
the mass density, the total velocity, the perturbed velocity and 
magnetic field (see Section [2TJ . We normalise the time with the 
orbital timescale so that Q = 1 in these time units. In addition, 
we use a linear cooling function A(T) = aT + b with a and b 
constant parameters and we write: 

^)- d( " 0) -A((e)). (12) 

We offset the mechanical energy m by its value mo (constant in 
time) on the equilibrium velocity profile. If we further neglect the 
potential energy in the mechanical energy, then 

(m- mo) ~ {\pu 2 + -!-b 2 }. (13) 



d(e) 
At 
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Hence both ( W xy ) and ( m — mo ) are quadratic in the perturba- 
tion amplitude. The comparison of equation (1121 with equation 
{7} leads us to use 



f =(e> 

as the reduced variable for temperature and 

1 



8tt 



(14) 



(15) 



as the reduced variable for the turbulent amplitude. 



3 RESULTS FOR CUBIC DOMAINS 

We select four initial conditions in each of the four quadrants that 
produce (f , y) simulation trajectories qualitatively similar to the 
(T, y) trajectories of the reduced model (thick coloured lines in 
Fig. 11. 



3.1 Parameters 

In order to make the comparison more quantitative, we need to 
specify the parameters A, n and W . We use the estimated posi- 
tions of the points T and 'A' in the phase space (f , y) in order 



Figure 1. Phase diagram of the reduced model equations (J9j 
and ITOt with parameters A = 0.1, n = 2 and W = 1.2. The 
separatrices of point T are shown as dotted and dashed lines. 



to constrain the set of parameters that yields the best agreement. 
Indeed, the coordinates (Tj, j/j) for the critical point T completely 
specify the parameters A and W for a given value of n: 



and 



W ■■ 



API) 



(16) 



(17) 



The point T is located at the crossing point between the 
separatrices, which are the asymptotes of the trajectories at early 
and late times. We estimate the position of the point T for the 
simulations from the point at which the red and green trajectories 
diverge (see Fig. |2(aJ| . 

In fact, the equations above also hold with A indices : the 
active saturation point can also be used to fit the parameters. We 
find that A = 0.1, n = 2 and W = 1.2 are values which satisfy 
both the conditions for points T and 'A'. 



3.2 Phase diagram comparison 



We plot the simulation trajectories in Fig. s [2(a)] and [2(b) | These 
should be compared with the trajectories in the first Fig. 1. As 
is plain, the qualitative agreement between the two methods is 
excellent. Particularly remarkable is the quantitative agreement 
near the critical point. This demonstrates that the relative values 
of thermal and turbulent time scales are accurately reproduced 
by the reduced model (the phase space trajectories are governed 
by the ratio between dy/dt and dT/di)). Fig. |2(b"J| also shows 
that the reduced model can describe the behaviour of the 3D 
system in a consistent way over a great range of magnitudes for 
y. The small value for the parameter A actually helps the linear 
approximation to hold until larger amplitudes for y. The fastest 
linear modes are then likely to grow to large amplitudes relative 
to the other mode s, which makes th i s syst em suitable to parasiti c 
analyses such as in iGoodman fc Xul < 19941) or lLatter et alJteOOSl) . 

Our reduced model is in effect a 2 variable projection of a 
system with many more degrees of freedom. But surprisingly few 
details show the limits of the reduced model: only some trajecto- 
ries exhibit self-crossing behaviour, for example when the sytem 
wanders near the saturation point 'A'. 
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Phase space (wpointA) 




„ 4 (to point Q) 6 , 

Temperature (T) 



(a) Close up on the point T. 



Phase space 
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(b) Global view of the phase diagram. 



1000.0 



Figure 2. Comparison with simulations of cubic domains. The 
separatrices of Fig. 1 are overlayed with five test trajectories (T, y) 
in our simulations. The dotted green trajectory is the case 'Ct' 
simulation. 



3.3 The route towards saturation 

Here we focus on a typical cubic domain simulation, namely, the 
case 'Ct' simulation (see appendix IA5I for the nomenclature: case 
'Ci' is isothermal, case 'Cc' has cooling, case 'Ct' has cooling and a 
temperature-dependent resistivity). But these comments are also 
valid for other cubic domain cases 'Ci' and 'Cc'. From the ini- 
tial conditions, the fastest growing mode is a vertical 2-channel 
flow which soon dominates. The mode grows in amplitude on a 
time scale no shorter than the orbital time scale. Meanwhile, the 
Alfven speed increases and its crossing time through the box gets 
shorter and shorter. The total pressure balance across the simula- 
tion box is hence better and better realised. Therefore the growth 
of the 2-channel proceeds with a nearly uniform total pressure. 
The location of the two counter flowing channels corresponds to 
the nulls of magnetic pressure. Consequently, the thermal pres- 
sure must be maximal within the channels, and that pressure 
must increase as the amplitude of the mode rises. The result is 
a sharp increase in density at the location of the channels which 
get thinner and thinner due to mass conservation. 

On either side of these channels, the magnetic field is almost 
uniform with same direction and opposite sign. The channels are 
also the current sheets where these two ordered fields reconnect. 
Ohmic heating is thus maximal at the location of these channels, 



which experience the highest temperature in the computational 
domain for non isothermal simulations (case 'c' or 't'). 

Fig. |3(b)| illustrates the simulation 'Ct' during this phase of 
channel growth. The surprising sharpness of the channel feature 
is due both to the cube geometry of the box and to the linear 
filtering efficiency of case 't'. Indeed, when the resistivity increases 
(for lower temperatures) towards the marginal stability, fewer and 
fewer modes are available until only one mode is unstable. Strong 
channel features are therefore likely to be a characteristic of flows 
with a temperature close to the critical T point. Also, the cubic 
geometry of t he box does not al l ow fo r unstable parasitic modes as 
computed bv lGoodman fe Xul jl994) to fit in the computational 
domain, so that the growing channel m ode remains unhind ered 
until very high amplitudes. As detailed in lLatter et all d2009l ') (for 
isothermal gases), it is then the shrinking of the channel width 
which decreases the wavelength of unstable parasites so that they 
finally fit in the box and start attacking the channel flow. 

The growing process of the channel flow can also be viewed 
as the progressive spatial (as opposed to temporal as in the phase 
portrait) separation of the system into two phases: a high tem- 
perature region of dense and fast moving fluid with a disordered 
magnetic field (within the channels), and a lower temperature re- 
gion of dilute quiet fluid with a strong ordered magnetic field (in 
between the channels). This is of course reminiscent of the phase 
separation described by our reduced model into an active and a 
quiet region. However, the two states are now coupled. For ex- 
ample, the fluid in the final quieter re gion does transpor t some 
angular momentum, as in the work of iFleming fc Stone! (|2003h 
where a dead zone is activated by its neighbouring active zone. It 
also has a temperature significantly higher than the quiet state 
temperature To due both to thermal diffusion and to the dissipa- 
tion of the magnetic field. The partitioning of the medium into 
a dense and dilute phase, while the total pressure is kept uni- 
form, resembles the classical Field thermal instability: high and 
low density regions appear at uniform thermal pressure. We recall 
that lBalbus fc Lesaffre] (|2008h interpreted the stability around the 
point T thanks to a criterion similar to the Field criterion. The 
present study provides further evidence for this view. 



3.4 The saturated state 

Even after the onset of the parasitic modes, both channels remain 
relatively well defined: we find that cubic domains lead to perma- 
nent channel flows (as witnessed in lLatter et al.|[2009l) . The par- 
asites soon saturate, and the two channels remain flowing with 
travelling wave-like structures (as in saturated kink modes) or 
occasional plasmoids (saturated pinch modes) as illustrated on 
Fig. |3(c)| These structures are reminiscent of the saturated kink 
and pinch modes th at feed on magnetised jets as observed in 
Bisk amp et al.l fll99d ). The parasites' perturbation of one channel 
are sometimes strong enough to diffuse and meet the parasites of 
the other channel, at which point some stronger coupling occurs, 
the channel flow slows and becomes more disordered, although 
two blocks with opposite velocities can still readily be identified 
as seen on Fig. |3(d)| These events correspond to sudden dips in 
the evolution of the total energy contained in the computational 
domain. After a short while, two well defined channels reform, 
with their saturated parasites. 

Simulations with shorter cooling time scales have less pres- 
sure support within the channel flows. They show thinner chan- 
nels which break down sooner and yield lower saturation levels. 



4 OTHER GEOMETRIES 

We have thus far limited the discussion to cubic domains, as these 
seem particularly fit for our reduced model. We briefly report here 
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(a) Initial conditions (after first 
transient evolution) 




(b) The growing channel 




(c) The steady state 




(d) At a trough 



Figure 3. Several density slices at various times of the simula- 
tion in case 'Ct'. The evolution of the total kinetic plus magnetic 
energy (y 2 ) in the perturbation is also given in the top right hand 
corner of each figure. 

some of the potentially interesting behaviour we encountered in 
other geometries. 

4.1 Bar results 

Many MRI simulations have been performed in compu- 
tational dom a ins el ongated in t h e az i muthal direction 
dHawley et alj Il995l: iFleming et al.l l2000l: ISano fc Inutsukal 
200ll; iFromane fc Papaloizo ui 120071 ; iFromang et al.l l2007t 
Lcsur fe Longarettil l2007h . Here we focus on the results we 
obtained in a similar geometry (our case 'B' with L x = 1 and 

Ly = 4). 

The bar geometry still does not allow for any 




Figure 4. Density slices for case 'BtO' on its way near the tip of 
a peak in total energy. 



iGoodman fc Xul lll994h parasites to fit in the computational 
domain for the channel flow of maximum growth. Indeed, the 
parasites' wave vectors are restricted to a sector of angle less than 
7r/2 encompassing the d irection of the channel (Goodma n fc Xul 
1 1994 lLatter et afll200Sh . The direction of the dominant channel 
flow is close to 7r/4 with respect to the direction axes. Hence the 
x-direction wavelength of a parasite needs to be bigger than the 
vertical extent of the channel in order to grow. As the channel 
gets thinner, however, compressible parasites ma y fit in the 
comp utational domain and attack the structure jLatter et al.l 

boosh . 

The chief difference between the cube and the bar is that in 
the latter there is enough room in the y-direction for the para- 
sites to divert one channel into the path of the other, an event 
that leads to a dramatic collapse of the whole channel structure. 
Non-linear interactions quickly render the flow isotropic while re- 
connection brings the perturbation back to linear amplitudes. De- 
pending on the overall resistivity (hence on the overall tempera- 
ture), a new channel mode can or cannot grow. Depending on its 
thermal inertia, the system either enters a new cycle of eruption 
or decays. We can get immediate decay (case 'Bt'), a few mild 
eruptions before decay (case 'Bt' with less initial seed for para- 
sites), infinite sequence of mild eruptions (case 'Bt' with a longer 
cooling timescale), infinite sequence of moderate eruptions (case 
'Bi' and 'Be') or infinite sequence of strong eruptions (case 'BtO', 
see Fig. [4}. In case 'BtO', the resistivity function is designed such 
that even for the radiative balance temperature there is still one 
unstable vertical MRI mode. In the next section, we devise a new 
reduced model which accounts for and sheds light on this wide 
variety of behaviour. 



4.2 Eruptive reduced model 

The previous discussion emphasises the role of the parasitic modes 
for the main growing mode which usually is a 2-channel mode. 
These paras i tes w ere described in IGoodman fc X u| lll994f) and 
lLatter et all 1120091) as a horizontal perturbation on top of a ver- 
tical main mode of evolution. Here wc separate each field variable 
X into its horizontally averaged contribution ( X ) z (which only 
depends on the altitude z) and its remainder SX = X — (X) z . 
The component ( X ) z naturally traces the main vertical 'channel' 
when one is dominating. And the component 8X can be thought 
of the 'perturbations' on top of this average profile. 

We now look at the separate evolution of the energies con- 
tained in each of the components. In particular, we monitor the 
magnetic energy in the 'channel' component 

y 2 c = ^((b) z -(b) z ) (18) 
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Thermal, 'channel' and 'perturbations' energies 
1000.000 




5 10 is 

Time (orbits) 




Figure 5. Evolution of variables T (red, solid), y% (blue, dashed) 
and j/p (green, dotted) (see text for definitions) in the case 'BtO' 
simulation. 



Figure 6. Evolution of variables T, and y p for the reduced 
equations p0 ) -p2 )l with A = 0.3, W = 3/2. The initial values 
are T = 1, y c = 1 and y p = 0.1. 




and in the 'perturbations' component 

Vl = ^(b 2 -(b) z .(b) z ) (19) 

which are displayed on Fig. [5] along with the quantity T = (e). 
We set the beginning of an eruptive cycle at a minimum of the 
'channel' energy, where a dominant channel mode is just starting 
to emerge and grow. At this point 'perturbations' are still decay- 
ing roughly at an orbital rate and they dominate the total en- 
ergy: the medium is vertically homogeneous in a statistical sense. 
The perturbations start to grow only after the channel amplitude 
has reached a certain threshold around 1 in our units. Then on 
the growth rate of the perturbations keeps increasing, in agree- 
ment with the idea tha t the growth of parasites is proportional t o 
the channel amplitude llGoodman fc Xjl994l:lLatter et al.l2009l) . 
When the perturbations reach an amplitude comparable to the 
channel, the channel collapses and the medium is perturbation 
dominated. At this point, the system is in a state of decaying 
turbulence superimposed with spiral density waves, until the am- 
plitudes become linear again, and a new channel mode emerges. 

The description above naturally leads to a simple set of re- 
duced equations where our former variable y is now detailed into 
two components y c and y p : 

^ = a(T)y c - Ay 2 c y p , (20) 
at 

^ = (Vc - l)»p (21) 

and 

^ = Wy 2 ~^-A(T). (22) 
at at 

The former non-linear saturation term is now replaced by a cou- 
pling term between the channel and the perturbations. The square 
exponent for y c in Ay^y p accounts for the observation that the 
breakup of a channel arises when two parasites interact. This sim- 
ple reduced model is able to reproduce qualitatively each of the 
behaviours experienced in the simulations of the previous section, 
provided it is started with the right set of initial conditions for the 
channel and perturbation component (see Fig. [6] for case 'BtO'). 

The main missing piece in this reduced model is how the 
decaying part of the cycle determines the starting amplitudes for 
the next cycle. The dice are rolled again during each phase of 
decaying turbulence. 



Figure 7. Density slices for the case 'Sc' simulation at a point 
where a channel-like structure stands out. 

4.3 Slab results 

We also tested geometries with a slab aspect ratio L x = L y = 4 
while L z = 1. In this case, iGoodman fc Xul lll994h parasites can 
fit in the computational domain, but the amplitude of the result- 
ing channel s remain small so that a parasite analysis cannot be 
applied (see lLatter et aLll2009T ). Case 'St' decays from the begin- 
ning, as for case 'Bt', but case 'StO' (with rji = 0.007) is similar 
to case 'Sc': intermittently, at the peak of the magnetic activity, a 
kinky channel-like structure (see Fi g.O stands out. It washes out 
very quickly. This agrees with what Bodo et al. (2008) observe in 
their aspect ratio study. 

The flow cannot be described as accurately with a simple 
model, as many modes of order 1 amplitude are present, all inter- 
acting with one another. However, the reduced model of equations 
||9} and 1101 may still be used for qualitative comparison, but with 
a higher value for the constant A. 



5 SUMMARY AND PROSPECTS. 

We have performed Cartesian shearing-box simulations in condi- 
tions near marginal stability, which exist at the border between 
dead and active zones in accretion disks. The linear relative fil- 
tering at low amplitudes helps select prominent channel flows. 
This makes the system suitable for a description with simple re- 
duced models which are able to account for general properties of 
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the simulations such as internal energy, mechanical energy and 
momentum transport. 

We investigated the influence of the thermal properties of 
the gas on the outcome of our simulations. The aspect ratio of 
our computational domain also has a very strong impact on their 
outcome. This is very likely due to the non- linear dynamics of 
interactions between various MRI modes (see lLatter et aTll2009l. 
for details). 

Stratified and global simulations of disks are now required 
to determine which geometrical configuration is selected in real 
disks. Also, the presence of a mean vertical magnetic field in 
our simulations sets a preferred scale for the channel: simulations 
without zero net mean field still need to be investigated. 

Last, we implied a dependence of the resistivity on the tem- 
perature which originally is mainly due to the electron fraction. 
The time scales for the recombination of electrons can be large 
and time-dependent chemistry has then to be accounted for. This 
eventually has to be i ncluded in order to link with the work of 
lllgner fc Ne lson (2008), for example. 



control is made non-uniform accordingly. Similarly, note that vis- 
cosity is not assumed uniform in the present implementation of 
the viscous term, although we use it as a constant in the present 
study. 

The diffusion step is done at the end of the standard ZEUS3D 
step and assumes the density is kept fixed. The cooling term fol- 
lows and also assumes the density does not change. The dissi- 
pative source heat terms e<j (viscous and Ohmic) are computed 
along with the resistive and viscous terms and added as a con- 
stant to the cooling function At = A + e^. We adopt a linear 
cooling function A(T) = aT + b which allows us to compute the 
isochore evolution of the internal energy analytically. Time-step 
controls for the thermal diffusion and cooling are based on the dif- 
fusion time across a pixel Ax 2 /x (where Ax is the size of a pixel) 
and the local cooling time-scale ^ \ These controls are never 
reached in practise (the limiting time is generally the Courant- 
Levy-Friedrich condition in the hottest pixel). 

The version of the code used in this paper is available for 
download on the MAGNET website [http: / /magnet .ens.fr 
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APPENDIX A: NUMERICAL SETUP 
Al Numerical method 



Our c ode is basically the ZEUS3D code llStone fc Norman 
[l992a|2) with the shearing-box implemented as in lHawIeT"et al 



1995), and we started from the version used by S. Fromang 



Fromang &: Papa loizou 200]}) but amended it to better suit the 



problem at hand. We slightly improved the parallelisation scheme 
to increase its efficiency and accuracy. 

We reset the mean magnetic field and the mean momen- 
tum t o their initial values every 150 time steps, as in ISilvers! 
(2008). This prevents truncation errors for the mean momentum 
and magnetic fields from wandering too far from zero. 

We split the momentum transport step for the Euler equa- 
tion. We first transport the perturbed velocity u. This increases 
the accuracy of the Van-Leer slopes. Then we combine the ad- 
vection of the mean shear with the Coriolis a nd tidal forces in 
an ex act analytic rotation as in Section 3.3 of Grcsscl & Zieglerl 
||2007| ). This increases the accuracy of the treatment of epicylic 
motions. 

We do not use the artificial viscosity of ZEUS3D. We rather 
implement a physical viscous term which takes the same form as 
in equation J2t . We rely on this term to degrade kinetic energy 
into heat and ensure total energy conservation (see below in Sec- 
tion lA5l l. We rewrite the resistive term in order to be able to use a 
non uniform resistivity r\. The resistive criterion for the time-step 



A2 Units 

We use the orbital timescale I/O as our unit time. We use the 
vertical extent L z of the computational domain as our unit length. 
We use the mass initially contained in a cube of size L z as our unit 
mass, so that the initial average density in the box is ( p ) = 1 (the 
brackets denote volume averages over the full box). In addition, 
most of our runs begin with (p) = 1 in our units. The initial 
pressure scale height is therefore H = 1/^/7 — 0.77L Z : the unit 
length is of the order of the initial pressure scale height. The scale 
height later varies as the average temperature in the box changes. 



A3 Thermal properties 

In the isothermal case, the pressure is simply related to the den- 
sity by p = CqP where cq is the initial sound speed in the non- 
isothermal runs. For the non-isothermal cases, the cooling func- 
tion A is a function of the temperature which we take to be linear 

A(T) = aT + b. 

In effect, this can be understood as the first order Taylor expan- 
sion of a more realistic cooling function around an intermediate 
temperature between the active and quiet states mentioned in the 
introduction. 

The steady-state temperature To results from the balance 
between viscous heating and radiative cooling. It may thus be 
written as 



T = (u 



q 2 Q 2 



b)/a 



q 2 U 2 
4a 



(Al) 



which is slightly higher than the radiative equilibrium temper- 
ature T r = —b/a due to viscous heating. We fix the radiative 
equilibrium temperature to T r = 1/2 and hence adopt 



1, 



A(T) = a(T--). 



(A2) 



The cooling time scale then depends only on the parame- 
ter a, which we set to a = 1 in most of our runs. Larger values 
for a yield similar results as in the isothermal case. Smaller val- 
ues mean a longer cooling time 1/a and consequently a statisti- 
cally stationary thermal state is reached after a longer integration 
time. Small values for a also yield higher equilibrium tempera- 
tures which constrain the time-step to much lower values due to 
the Courant-Friedrich-Levy (CFL) condition. 
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Figure Al. MRI dispersion relation for the vertical wave num- 
bers in the conditions of the main simulations of this paper. Vis- 
cosity and resistivity are included. We display the real part of 
the fastest growing eigenvalue in units of the orbital time scale. 
The actual dispersion relation can be found in Lesaffre & Balbus 
(2007) , equation (16). The dotted line is for no ideal MHD. The 
solid line is for case 'c' (tj = v = 2 X 10 — 3 ). The dashed line is 
for case 't' with T = 0.5 (radiative equilibrium) and the dash- 
dotted line is for T = 2/3 (as in the initial conditions, which have 
(p) = 1). The vertical line indicates the vertical wavenumber of 
the box. The wavenumbers resolved by our simulations are right 
of this line, up to k = 647r for our standard resolution. 



A4 Initial conditions 

The mean vertical magnetic field is set to the value Bo = 
^/2/4007r which corresponds to a plasma f} value of around 10 2 
when the average pressure is ( p ) = 1 . For that value of the mag- 
netic field, the wavelength of the MRI mode of maximum growth 
rate just fits in the vertical extent of the box (see Fig. IAl| . 

Because the system supports two stationary states accessed 
by different initial conditions, we must then be careful how we 
characterise these. The initial configuration consists of a pertur- 
bation added to the stationary flow described above. 

Ra ndom perturbations with low initial amplitude such as 
used bv lHawlev et al] <|l995h lead to a long initial transient phase 
during which MRI modes grow and other components of the per- 
turbation field decay, until the fastest growing MRI modes stand 
out. During this linear filtering phase, the temperature as well as 
the magnitude of the perturbation may change significantly. 

In order to shorten this phase, we prepare the system in a 
state that is already linearly filtered to some extent. We choose to 
build the initial perturbation with a combination of some of the 
fastest growing modes and add a random seed. That mixture of a 
random perturbation and a selection of growing modes consider- 
ably shortens the initial transient phase. With this setup, we are 
indeed able to eliminate the initial strong first channel flow which 
appears in the time evolution of MRI shearing-box simulations. 

We set the initial perturbation in the Fourier space of a larger 
domain and with a small scale cut off at intermediate scales. We 
then select only those modes which fit in the box. That way the 
same perturbation is used for runs with different geometries, dif- 
ferent resolutions or different ways of sharing the domain between 
the processors units. To work in the Fourier space also allows to 
properly clean the divergence of the initial magnetic field. 

Finally, we normalise the initial perturbation such that the 
maximum radial velocity is sonic (in fact, equal to 1 in our units). 
That high initial amplitude maximises our chances to first pick 
up the active state. 



A5 Parameter study 

The main parameter we varied was the aspect ratio of the com- 
putational domain. The vertical size of the box L z was kept fixed 
to L z = 1. We set its radial size L x and azimuthal size L y to 
L x = L y = 1 for the cube (first letter 'C'), to L x = 1 and L y = 4 
for a more classical azimuthal bar (first letter 'B') aspect ratio, 
and to L x = L y = 4 for a slab (first letter 'S') aspect ratio. 

We also changed the microphysical properties of the gas. We 
tested isothermal runs (second letter 'i'), runs with constant re- 
sistivity and cooling (second letter 'c') and runs with a varying 
resistivity and cooling (second letter 't'). Some runs were per- 
formed with a lower resistivity and have a trailing zero in their 
coding name (simulation 'BtO', for example, is a bar simulation 
with a varying resistivity lower than that of simulation 'Bt'). 



A6 Resolution 

We us e an isotropic resolution as advised by iLesaffre fc Balbus! 
(2007). Our pixels are cubes of side length Ax = 1/64 for our 
standard resolution. We halved that value for the case 't' of each 
of the geometries for a resolution test on shorter integration times. 
The size of our simulations was then, in terms of zone numbers: 
64 3 for the cubes, 64 X 256 X 64 for the bars and 256 X 256 X 64 for 
the slabs (and eight times these numbers for the higher resolu- 
tion t ests) . Our value for the Courant number llLesaffre fc Balbus! 
|200?t ) ranges from 0.2 for cases 'B' and 'S' to 0.5 for case 'C. A 
low Courant number is necessary especially for case 't' where 
large resistivity values can be obtained because our scheme uses 
an explicit resistive term. We aimed at fulfilling an integration 
time of a hundred of orbits. However, that time span could not 
be reached in case 'C Case 'BtO' is run for up to 330 orbits. 
In general, our simulations are run for time-step numbers of the 
order of a million up to six millions. 



A7 Dissipation and diffusion coefficients 

Total energy is best conserved when the physical dissipation is 
much h igher than the numerica l dissipation. We use equation 
(50) of ILesaffre fc Balbus! [120071 ) with Ai = 1/64, /3 = 100, 
C = 0.5 and k = 64-7T to get an upper bound estimate of the total 
numerical dissipation coefficient r/-^ < 5x 10~ 3 . Accord- 

ingly, we choose our minimal physical dissipation coefficients as 
rjo = v = 2 X 10 -3 . Our standard runs are hence only marginally 
resolved in places where strong gradients occur or where the mag- 
netic field becomes very large. Our control runs with twice the 
resolution show that this does not impact our results except in 
case 'C, where we show only the results at a resolution of 128 3 . 
According to the same line of thought, we include thermal dif- 
fusion with \ = 2 X 10 — 3 in cases 'c' and 't' in order to avoid 
spurious thermal diffusion effects. 

The resistivity in cases 'i' and 'c' is simply set to r\ = r/rj- In 
case 't', we adopted a resistivity of the form 

r,(T)= Vo + m T-^ 2 (A3) 

to mimic the power law dependence of the resistivity in a fully 
ionised plasma. The purpose was to test the validity of our sim- 
ple dynamical system: both an active and a quiet state should be 
within reach of the system. Hence, the coefficient 771 is chosen rea- 
sonably small (to avoid putting too much constraint on the time 
step due to high values of resistivity), while forcing all the MRI 
modes that fit in the computational to decay when the tempera- 
ture is T = 0.5. We use the values 771 = 1.2 X 10 — 2 for standard 
runs and 771 = 7 X 10 -3 for lower resistivity runs (case 'BtO', for 
example). The MRI dispersion relations for various realisations 
of the dissipation coefficients are plotted in Fig. IA1I 
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